using TrixiParticles
particle_spacing = 0.05
domain_length = 1.0

function pressure_function(pos)
    t = 0
    U = 1.0
    b = -8pi^2 / 10
    x = pos[1]
    y = pos[2]

    return -U^2 * exp(2 * b * t) * (cos(4pi * x) + cos(4pi * y)) / 4
end

function velocity_function(pos)
    t = 0
    U = 1.0
    b = -8pi^2 / 10
    x = pos[1]
    y = pos[2]

    vel = U * exp(b * t) *
          [-cos(2pi * x) * sin(2pi * y), sin(2pi * x) * cos(2pi * y)]

    return SVector{2}(vel)
end

n_particles_xy = round(Int, domain_length / particle_spacing)

domain_fluid = RectangularShape(particle_spacing, (2, 1) .* n_particles_xy,
                                (0.0, 0.0), density=1000.0,
                                pressure=pressure_function,
                                velocity=velocity_function)

smoothing_length = 1.5 * particle_spacing
smoothing_kernel = WendlandC2Kernel{ndims(domain_fluid)}()
fluid_system = EntropicallyDampedSPHSystem(domain_fluid, smoothing_kernel,
                                           smoothing_length, 1.0)

fluid_system.cache.density .= domain_fluid.density

plane_out = ([2.0, 0.0], [2.0, domain_length])

outflow = BoundaryZone(; plane=plane_out, boundary_type=OutFlow(),
                       plane_normal=[-1.0, 0.0],
                       open_boundary_layers=10, density=1000.0, particle_spacing)
open_boundary_out = OpenBoundarySPHSystem(outflow; fluid_system,
                                          boundary_model=BoundaryModelTafuni(),
                                          buffer_size=0)

semi = Semidiscretization(fluid_system, open_boundary_out)
TrixiParticles.initialize_neighborhood_searches!(semi)

v_open_boundary = zero(outflow.initial_condition.velocity)
v_fluid = vcat(domain_fluid.velocity, domain_fluid.pressure')

TrixiParticles.set_zero!(open_boundary_out.pressure)

TrixiParticles.extrapolate_values!(open_boundary_out, v_open_boundary, v_fluid,
                                   outflow.initial_condition.coordinates,
                                   domain_fluid.coordinates, semi, 0.0;
                                   prescribed_pressure=false,
                                   prescribed_velocity=false)

trixi2vtk(fluid_system.initial_condition, filename="fluid",
          v=domain_fluid.velocity, p=domain_fluid.pressure)

trixi2vtk(open_boundary_out.initial_condition, filename="open_boundary_out",
          v=v_open_boundary, p=open_boundary_out.pressure)

plane_in = ([0.0, 0.0], [0.0, domain_length])

inflow = BoundaryZone(; plane=plane_in, boundary_type=InFlow(),
                       plane_normal=[1.0, 0.0],
                       open_boundary_layers=10, density=1000.0, particle_spacing)
open_boundary_in = OpenBoundarySPHSystem(inflow; fluid_system,
                                          boundary_model=BoundaryModelTafuni(),
                                          buffer_size=0)

semi = Semidiscretization(fluid_system, open_boundary_in)
TrixiParticles.initialize_neighborhood_searches!(semi)

v_open_boundary = zero(inflow.initial_condition.velocity)
v_fluid = vcat(domain_fluid.velocity, domain_fluid.pressure')

TrixiParticles.set_zero!(open_boundary_in.pressure)

TrixiParticles.extrapolate_values!(open_boundary_in, v_open_boundary, v_fluid,
                                   inflow.initial_condition.coordinates,
                                   domain_fluid.coordinates, semi, 0.0;
                                   prescribed_pressure=false,
                                   prescribed_velocity=false)

trixi2vtk(fluid_system.initial_condition, filename="fluid",
          v=domain_fluid.velocity, p=domain_fluid.pressure)

trixi2vtk(open_boundary_in.initial_condition, filename="open_boundary_in",
          v=v_open_boundary, p=open_boundary_in.pressure)